Rclub_homework_Dec_12_2016
Stacey Harmer
December 12, 2016
Problem 1 : Tomato
Remember the tomato data set generated by Pepe; we first looked at this when we were working on Chapter 9. 35 accessions for seven species were grown in sun and shade.
Assess whether there is evidence for total length (“totleng”) response to shade and whether this response varies by species. Consider whether including accession (“acs”), and/or shelf (“shelf”) using adaptive priors improves the model fit.
Bonus: would it be better to consider shade by accession interactions instead of shade x species?
setwd("~/R /R_club/R_Club_2016/Stacey.Harmer_RethinkingHomework/Assigment_2016-12-12")
tom <- read.csv("TomatoR2CSHL.csv")
head(tom)## shelf flat col row acs trt days date hyp int1 int2 int3 int4
## 1 Z 1 B 1 LA2580 H 28 5/5/08 19.46 2.37 1.59 1.87 0.51
## 2 Z 1 C 1 LA1305 H 28 5/5/08 31.28 3.34 0.01 9.19 1.62
## 3 Z 1 D 1 LA1973 H 28 5/5/08 56.65 8.43 2.39 6.70 3.69
## 4 Z 1 E 1 LA2748 H 28 5/5/08 35.18 0.56 0.00 1.60 0.61
## 5 Z 1 F 1 LA2931 H 28 5/5/08 35.32 0.82 0.02 1.49 0.46
## 6 Z 1 G 1 LA1317 H 28 5/5/08 28.74 1.07 6.69 5.72 4.76
## intleng totleng petleng leafleng leafwid leafnum ndvi lat lon
## 1 6.34 25.80 15.78 30.53 34.44 5 111 -9.5167 -78.0083
## 2 14.16 45.44 12.36 22.93 13.99 4 120 -13.3833 -75.3583
## 3 21.21 77.86 13.05 46.71 43.78 5 110 -16.2333 -71.7000
## 4 2.77 37.95 8.08 26.82 33.28 5 105 -20.4833 -69.9833
## 5 2.79 38.11 7.68 22.40 23.61 5 106 -20.9167 -69.0667
## 6 18.24 46.98 23.66 42.35 42.35 5 132 -13.4167 -73.8417
## alt species who
## 1 740 S. pennellii Dan
## 2 3360 S. peruvianum Dan
## 3 2585 S. peruvianum Dan
## 4 1020 S. chilense Dan
## 5 2460 S. chilense Dan
## 6 2000 S. chmielewskii Danstr(tom)## 'data.frame': 1008 obs. of 25 variables:
## $ shelf : Factor w/ 6 levels "U","V","W","X",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ flat : int 1 1 1 1 1 1 1 1 1 1 ...
## $ col : Factor w/ 9 levels "A","B","C","D",..: 2 3 4 5 6 7 9 2 3 4 ...
## $ row : int 1 1 1 1 1 1 1 2 2 2 ...
## $ acs : Factor w/ 36 levels "LA1028","LA1305",..: 23 2 19 26 30 5 14 32 7 34 ...
## $ trt : Factor w/ 2 levels "H","L": 1 1 1 1 1 1 1 1 1 1 ...
## $ days : int 28 28 28 28 28 28 28 28 28 28 ...
## $ date : Factor w/ 2 levels "5/5/08","5/6/08": 1 1 1 1 1 1 1 1 1 1 ...
## $ hyp : num 19.5 31.3 56.6 35.2 35.3 ...
## $ int1 : num 2.37 3.34 8.43 0.56 0.82 1.07 2.85 2.08 5.43 4.08 ...
## $ int2 : num 1.59 0.01 2.39 0 0.02 6.69 0.41 0.53 0.81 3.26 ...
## $ int3 : num 1.87 9.19 6.7 1.6 1.49 5.72 3.79 1.9 3.63 3.49 ...
## $ int4 : num 0.51 1.62 3.69 0.61 0.46 4.76 3.25 NA 2.66 3.02 ...
## $ intleng : num 6.34 14.16 21.21 2.77 2.79 ...
## $ totleng : num 25.8 45.4 77.9 38 38.1 ...
## $ petleng : num 15.78 12.36 13.05 8.08 7.68 ...
## $ leafleng: num 30.5 22.9 46.7 26.8 22.4 ...
## $ leafwid : num 34.4 14 43.8 33.3 23.6 ...
## $ leafnum : int 5 4 5 5 5 5 5 4 5 5 ...
## $ ndvi : int 111 120 110 105 106 132 118 112 107 123 ...
## $ lat : num -9.52 -13.38 -16.23 -20.48 -20.92 ...
## $ lon : num -78 -75.4 -71.7 -70 -69.1 ...
## $ alt : int 740 3360 2585 1020 2460 2000 2920 480 75 3540 ...
## $ species : Factor w/ 5 levels "S. chilense",..: 4 5 5 1 1 2 3 4 5 5 ...
## $ who : Factor w/ 2 levels "Dan","Pepe": 1 1 1 1 1 1 1 1 1 1 ...summary(tom)## shelf flat col row acs
## U:161 Min. : 1.00 G :133 Min. :1.00 LA1954 : 40
## V:174 1st Qu.: 9.00 H :127 1st Qu.:2.00 LA2695 : 39
## W:178 Median :17.00 F :125 Median :3.00 LA1361 : 37
## X:174 Mean :17.89 C :117 Mean :2.56 LA2167 : 37
## Y:125 3rd Qu.:28.00 D :117 3rd Qu.:4.00 LA2773 : 37
## Z:196 Max. :36.00 E :107 Max. :4.00 LA1474 : 36
## (Other):282 (Other):782
## trt days date hyp int1
## H:495 Min. :28.00 5/5/08:716 Min. : 6.17 Min. : 0.00
## L:513 1st Qu.:28.00 5/6/08:292 1st Qu.:26.81 1st Qu.: 1.74
## Median :28.00 Median :32.02 Median : 3.59
## Mean :28.29 Mean :33.36 Mean : 4.71
## 3rd Qu.:29.00 3rd Qu.:38.56 3rd Qu.: 6.46
## Max. :29.00 Max. :74.60 Max. :39.01
## NA's :1
## int2 int3 int4 intleng
## Min. : 0.000 Min. : 0.010 Min. : 0.030 Min. : 0.000
## 1st Qu.: 1.060 1st Qu.: 2.975 1st Qu.: 2.163 1st Qu.: 9.637
## Median : 3.120 Median : 5.625 Median : 3.995 Median :17.255
## Mean : 4.287 Mean : 6.794 Mean : 5.102 Mean :20.340
## 3rd Qu.: 6.320 3rd Qu.: 9.367 3rd Qu.: 7.018 3rd Qu.:28.145
## Max. :28.980 Max. :27.760 Max. :23.280 Max. :92.420
## NA's :1 NA's :4 NA's :102
## totleng petleng leafleng leafwid
## Min. : 13.59 Min. : 1.53 Min. : 9.74 Min. : 8.29
## 1st Qu.: 39.25 1st Qu.:11.20 1st Qu.:27.43 1st Qu.:29.48
## Median : 50.98 Median :15.13 Median :34.59 Median :39.62
## Mean : 53.70 Mean :15.92 Mean :35.54 Mean :39.29
## 3rd Qu.: 64.76 3rd Qu.:20.48 3rd Qu.:42.98 3rd Qu.:47.75
## Max. :129.43 Max. :44.44 Max. :95.19 Max. :90.27
## NA's :2 NA's :1 NA's :1
## leafnum ndvi lat lon
## Min. :3.000 Min. :100 Min. :-25.400 Min. :-78.52
## 1st Qu.:5.000 1st Qu.:108 1st Qu.:-16.607 1st Qu.:-75.92
## Median :5.000 Median :115 Median :-14.152 Median :-73.63
## Mean :5.063 Mean :118 Mean :-14.490 Mean :-73.71
## 3rd Qu.:6.000 3rd Qu.:128 3rd Qu.:-12.450 3rd Qu.:-71.70
## Max. :8.000 Max. :137 Max. : -5.767 Max. :-68.07
## NA's :1
## alt species who
## Min. : 0 S. chilense :207 Dan :402
## 1st Qu.:1020 S. chmielewskii:226 Pepe:606
## Median :2240 S. habrochaites:226
## Mean :2035 S. pennellii :132
## 3rd Qu.:3110 S. peruvianum :217
## Max. :3540
## library(rethinking)## Loading required package: rstan## Loading required package: ggplot2## Loading required package: StanHeaders## rstan (Version 2.12.1, packaged: 2016-09-11 13:07:50 UTC, GitRev: 85f7a56811da)## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())## Loading required package: parallel## rethinking (Version 1.59)# get map2stan up and ready
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())Make dummy variable, so that sun = 0, shade = 1
tom$shade <- ifelse(tom$trt == "L", 1, 0)
head(tom, n=6)## shelf flat col row acs trt days date hyp int1 int2 int3 int4
## 1 Z 1 B 1 LA2580 H 28 5/5/08 19.46 2.37 1.59 1.87 0.51
## 2 Z 1 C 1 LA1305 H 28 5/5/08 31.28 3.34 0.01 9.19 1.62
## 3 Z 1 D 1 LA1973 H 28 5/5/08 56.65 8.43 2.39 6.70 3.69
## 4 Z 1 E 1 LA2748 H 28 5/5/08 35.18 0.56 0.00 1.60 0.61
## 5 Z 1 F 1 LA2931 H 28 5/5/08 35.32 0.82 0.02 1.49 0.46
## 6 Z 1 G 1 LA1317 H 28 5/5/08 28.74 1.07 6.69 5.72 4.76
## intleng totleng petleng leafleng leafwid leafnum ndvi lat lon
## 1 6.34 25.80 15.78 30.53 34.44 5 111 -9.5167 -78.0083
## 2 14.16 45.44 12.36 22.93 13.99 4 120 -13.3833 -75.3583
## 3 21.21 77.86 13.05 46.71 43.78 5 110 -16.2333 -71.7000
## 4 2.77 37.95 8.08 26.82 33.28 5 105 -20.4833 -69.9833
## 5 2.79 38.11 7.68 22.40 23.61 5 106 -20.9167 -69.0667
## 6 18.24 46.98 23.66 42.35 42.35 5 132 -13.4167 -73.8417
## alt species who shade
## 1 740 S. pennellii Dan 0
## 2 3360 S. peruvianum Dan 0
## 3 2585 S. peruvianum Dan 0
## 4 1020 S. chilense Dan 0
## 5 2460 S. chilense Dan 0
## 6 2000 S. chmielewskii Dan 0tail(tom, n=20)## shelf flat col row acs trt days date hyp int1 int2 int3 int4
## 989 X 36 A 2 LA2167 H 28 5/5/08 30.00 1.71 3.00 4.16 4.46
## 990 X 36 B 2 LA3788 H 28 5/5/08 19.70 0.62 0.74 0.90 0.43
## 991 X 36 C 2 LA1336 H 28 5/5/08 35.19 3.43 2.64 5.93 5.27
## 992 X 36 D 2 LA3795 H 28 5/5/08 26.17 2.06 1.99 0.87 5.65
## 993 X 36 E 2 LA2773 H 28 5/5/08 13.02 1.50 1.13 2.79 2.78
## 994 X 36 F 2 LA3788 H 28 5/5/08 23.39 5.47 7.44 9.95 4.74
## 995 X 36 H 2 LA1361 H 28 5/5/08 26.34 1.75 2.46 5.42 3.88
## 996 X 36 A 3 LA1969 H 28 5/5/08 26.36 0.16 0.61 0.30 0.13
## 997 X 36 C 3 LA1474 H 28 5/5/08 34.95 2.87 2.58 5.27 4.28
## 998 X 36 E 3 LA2884 H 28 5/5/08 26.32 1.35 0.44 3.46 2.11
## 999 X 36 F 3 LA1306 H 28 5/5/08 27.06 1.86 5.96 4.97 3.00
## 1000 X 36 G 3 LA2663 H 28 5/5/08 26.03 1.42 4.27 5.35 3.97
## 1001 X 36 H 3 LA1363 H 28 5/5/08 30.61 3.87 4.44 4.49 2.90
## 1002 X 36 I 3 LA2167 H 28 5/5/08 25.04 4.30 2.16 2.18 NA
## 1003 X 36 A 4 LA1474 H 28 5/5/08 29.38 0.12 0.44 0.44 0.21
## 1004 X 36 D 4 LA1969 H 28 5/5/08 37.00 1.01 0.72 0.70 1.20
## 1005 X 36 F 4 LA1316 H 28 5/5/08 26.59 2.00 5.17 6.06 3.01
## 1006 X 36 G 4 LA2695 H 28 5/5/08 32.11 1.26 2.29 3.69 4.55
## 1007 X 36 H 4 LA2695 H 28 5/5/08 29.00 3.06 7.37 6.14 4.12
## 1008 X 36 I 4 LA2204 H 28 5/5/08 27.62 2.03 4.99 3.75 NA
## intleng totleng petleng leafleng leafwid leafnum ndvi lat
## 989 13.33 43.33 12.42 49.90 51.83 6 124 -7.1800
## 990 2.69 22.39 8.03 23.08 24.92 4 112 -16.1383
## 991 17.27 52.46 16.86 33.21 43.20 6 107 -16.2094
## 992 10.57 36.74 NA 21.63 16.77 4 123 -10.1500
## 993 8.20 21.22 10.34 26.83 25.09 5 115 -18.2667
## 994 27.60 50.99 21.52 45.31 50.90 4 112 -16.1383
## 995 13.51 39.85 16.08 22.49 26.75 5 111 -9.5500
## 996 1.20 27.56 3.25 20.38 19.19 5 124 -17.5333
## 997 15.00 49.95 16.87 40.00 55.34 6 113 -16.6072
## 998 7.36 33.68 9.06 37.37 39.86 5 108 -23.6833
## 999 15.79 42.85 13.33 31.02 36.17 4 130 -12.9500
## 1000 15.01 41.04 16.12 40.15 48.18 5 137 -13.7333
## 1001 15.70 46.31 15.32 31.76 40.30 4 121 -10.1417
## 1002 8.64 33.68 13.94 16.96 22.22 5 124 -7.1800
## 1003 1.21 30.59 4.35 14.56 14.19 4 113 -16.6072
## 1004 3.63 40.63 6.42 18.67 23.17 5 124 -17.5333
## 1005 16.24 42.83 16.48 30.80 30.52 4 135 -13.3925
## 1006 11.79 43.90 23.77 42.93 50.53 5 128 -14.1525
## 1007 20.69 49.69 15.08 38.75 47.90 5 128 -14.1525
## 1008 10.77 38.39 10.25 22.97 26.65 4 128 -5.7667
## lon alt species who shade
## 989 -78.5200 1340 S. habrochaites Pepe 0
## 990 -73.6428 480 S. pennellii Pepe 0
## 991 -73.6222 75 S. peruvianum Pepe 0
## 992 -77.3500 3540 S. peruvianum Pepe 0
## 993 -69.5833 3260 S. chilense Pepe 0
## 994 -73.6428 480 S. pennellii Pepe 0
## 995 -77.8917 1340 S. habrochaites Pepe 0
## 996 -70.0333 3015 S. chilense Pepe 0
## 997 -72.6964 27 S. peruvianum Pepe 0
## 998 -68.0667 2390 S. chilense Pepe 0
## 999 -74.0167 3340 S. chmielewskii Pepe 0
## 1000 -71.9833 3110 S. chmielewskii Pepe 0
## 1001 -77.3833 3140 S. habrochaites Pepe 0
## 1002 -78.5200 1340 S. habrochaites Pepe 0
## 1003 -72.6964 27 S. peruvianum Pepe 0
## 1004 -70.0333 3015 S. chilense Pepe 0
## 1005 -73.9150 3100 S. chmielewskii Pepe 0
## 1006 -71.3886 3480 S. chmielewskii Pepe 0
## 1007 -71.3886 3480 S. chmielewskii Pepe 0
## 1008 -77.8667 2160 S. habrochaites Pepe 0I want to treat species as a cluster variable. Do I need to give each species a number? Try without doing that.
First write model only taking treatment and species into account After first attempt, learned that I need to trim this dataset. Keep accession (acs) and shelf and totleng and treatment and speices but ditch the rest
tom.trim <- tom[, c(1,5, 15, 24,26)]
head(tom.trim)## shelf acs totleng species shade
## 1 Z LA2580 25.80 S. pennellii 0
## 2 Z LA1305 45.44 S. peruvianum 0
## 3 Z LA1973 77.86 S. peruvianum 0
## 4 Z LA2748 37.95 S. chilense 0
## 5 Z LA2931 38.11 S. chilense 0
## 6 Z LA1317 46.98 S. chmielewskii 0str(tom.trim)## 'data.frame': 1008 obs. of 5 variables:
## $ shelf : Factor w/ 6 levels "U","V","W","X",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ acs : Factor w/ 36 levels "LA1028","LA1305",..: 23 2 19 26 30 5 14 32 7 34 ...
## $ totleng: num 25.8 45.4 77.9 38 38.1 ...
## $ species: Factor w/ 5 levels "S. chilense",..: 4 5 5 1 1 2 3 4 5 5 ...
## $ shade : num 0 0 0 0 0 0 0 0 0 0 ...tom.trim$sp_index <- coerce_index(tom.trim$species)mTom.1 <- map2stan(
alist(
totleng ~ dnorm(mu, sigma),
mu <- a_species[sp_index] + b_species[sp_index]*shade,
a_species[sp_index] ~ dnorm(0, 10),
b_species[sp_index] ~ dnorm(0,10),
sigma ~ dnorm(0, 10)), data=tom.trim, iter=5000, warmup=2000, chains =2)## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
## Elapsed Time: 1.75539 seconds (Warm-up)
## 1.36139 seconds (Sampling)
## 3.11678 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
## Elapsed Time: 1.55028 seconds (Warm-up)
## 1.34519 seconds (Sampling)
## 2.89547 seconds (Total)## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000142 seconds (Sampling)
## 0.000146 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]plot(mTom.1)
precis(mTom.1, depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_species[1] 36.88 1.40 34.66 39.10 4317 1
## a_species[2] 43.83 1.33 41.77 45.98 4190 1
## a_species[3] 44.34 1.30 42.26 46.42 4272 1
## a_species[4] 29.07 1.82 26.18 32.04 4082 1
## a_species[5] 48.20 1.36 46.02 50.33 3977 1
## b_species[1] 28.03 1.91 25.07 31.15 4164 1
## b_species[2] 17.10 1.83 14.10 19.86 4203 1
## b_species[3] 20.40 1.88 17.29 23.32 4192 1
## b_species[4] 26.62 2.42 22.74 30.41 4197 1
## b_species[5] 24.66 1.89 21.79 27.77 4522 1
## sigma 15.03 0.31 14.53 15.52 6000 1plot(precis(mTom.1, depth = 2))
Consider whether including accession (“acs”), and/or shelf (“shelf”) using adaptive priors improves the model fit.
# OK, include accession and then shelf; first make them an index
tom.trim$acs_index <- coerce_index(tom.trim$acs)
tom.trim$shelf_index <- coerce_index(tom.trim$shelf)
head(tom.trim)## shelf acs totleng species shade sp_index acs_index
## 1 Z LA2580 25.80 S. pennellii 0 4 23
## 2 Z LA1305 45.44 S. peruvianum 0 5 2
## 3 Z LA1973 77.86 S. peruvianum 0 5 19
## 4 Z LA2748 37.95 S. chilense 0 1 26
## 5 Z LA2931 38.11 S. chilense 0 1 30
## 6 Z LA1317 46.98 S. chmielewskii 0 2 5
## shelf_index
## 1 6
## 2 6
## 3 6
## 4 6
## 5 6
## 6 6# include accession, using adaptive prior
mTom.2 <- map2stan(
alist(
totleng ~ dnorm(mu, sigma),
mu <- a_species[sp_index] + a_acs[acs_index] + b_species[sp_index]*shade,
a_species[sp_index] ~ dnorm(0, 10),
b_species[sp_index] ~ dnorm(0,10),
a_acs[acs_index] ~ dnorm(0, sigma_acs),
sigma_acs ~ dnorm(0, 10),
sigma ~ dnorm(0, 10)), data=tom.trim, iter=5000, warmup=2000, chains =2)## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
## Elapsed Time: 6.35734 seconds (Warm-up)
## 5.93577 seconds (Sampling)
## 12.2931 seconds (Total)## The following numerical problems occured the indicated number of times after warmup on chain 1## count
## Exception thrown at line 21: normal_log: Scale parameter is 0, but must be > 0! 2## When a numerical problem occurs, the Hamiltonian proposal gets rejected.## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected## If the number in the 'count' column is small, do not ask about this message on stan-users.##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
## Elapsed Time: 5.7553 seconds (Warm-up)
## 4.98947 seconds (Sampling)
## 10.7448 seconds (Total)## The following numerical problems occured the indicated number of times after warmup on chain 2## count
## Exception thrown at line 21: normal_log: Scale parameter is 0, but must be > 0! 2## When a numerical problem occurs, the Hamiltonian proposal gets rejected.## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected## If the number in the 'count' column is small, do not ask about this message on stan-users.## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.000202 seconds (Sampling)
## 0.000205 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]plot(mTom.2)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
# neff not great for the a_acs values
precis(mTom.2, depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_species[1] 32.80 3.80 27.24 39.24 1146 1
## a_species[2] 39.05 3.93 33.09 45.30 1015 1
## a_species[3] 39.53 3.66 33.79 45.32 879 1
## a_species[4] 26.24 3.72 20.41 32.04 1214 1
## a_species[5] 42.67 3.83 37.18 48.86 1090 1
## b_species[1] 27.67 1.82 24.96 30.72 4889 1
## b_species[2] 16.83 1.74 14.07 19.56 4961 1
## b_species[3] 19.98 1.73 17.23 22.75 4424 1
## b_species[4] 26.98 2.26 23.70 30.86 4112 1
## b_species[5] 24.35 1.76 21.56 27.17 5099 1
## a_acs[1] -1.06 4.51 -7.83 6.32 1404 1
## a_acs[2] -6.38 4.24 -13.13 0.13 1412 1
## a_acs[3] 1.10 4.25 -5.96 7.24 1179 1
## a_acs[4] 7.99 4.35 0.78 14.34 1187 1
## a_acs[5] 18.75 4.40 11.49 25.19 1134 1
## a_acs[6] 7.29 4.36 0.63 14.21 1204 1
## a_acs[7] 14.95 4.34 8.18 21.63 1300 1
## a_acs[8] 10.98 4.08 4.88 17.55 1173 1
## a_acs[9] 8.03 4.11 1.50 14.44 982 1
## a_acs[10] 8.24 4.52 1.28 15.59 1648 1
## a_acs[11] 7.17 4.22 0.51 13.60 1269 1
## a_acs[12] -5.72 4.09 -12.36 0.49 1410 1
## a_acs[13] 0.17 4.18 -6.72 6.35 1508 1
## a_acs[14] 3.65 4.29 -2.98 10.65 1297 1
## a_acs[15] -0.65 5.43 -9.12 8.15 3129 1
## a_acs[16] 6.23 4.19 -0.80 12.32 1247 1
## a_acs[17] 9.35 4.20 2.61 15.76 1380 1
## a_acs[18] -8.07 4.08 -14.19 -1.46 1527 1
## a_acs[19] 20.27 4.46 13.20 27.16 1169 1
## a_acs[20] -1.86 3.99 -7.98 4.61 1197 1
## a_acs[21] 5.47 4.02 -0.46 12.28 1250 1
## a_acs[22] 10.57 5.15 2.18 18.52 1880 1
## a_acs[23] 2.00 4.23 -4.99 8.25 1587 1
## a_acs[24] 0.36 4.25 -5.87 7.44 1189 1
## a_acs[25] -0.46 4.21 -7.30 5.91 1174 1
## a_acs[26] 4.07 4.38 -2.33 11.49 1481 1
## a_acs[27] 6.38 4.15 -0.44 12.61 1334 1
## a_acs[28] -0.33 4.16 -6.99 6.19 1468 1
## a_acs[29] 5.56 4.31 -1.43 12.14 1447 1
## a_acs[30] 11.45 4.28 4.86 18.31 1354 1
## a_acs[31] 15.15 4.27 8.57 22.11 1213 1
## a_acs[32] 4.85 4.18 -1.96 11.16 1407 1
## a_acs[33] 0.94 4.11 -5.09 7.89 1501 1
## a_acs[34] -4.13 4.22 -10.99 2.14 1337 1
## a_acs[35] -1.17 4.31 -8.31 5.15 1418 1
## a_acs[36] -3.99 4.79 -11.94 3.36 2180 1
## sigma_acs 9.25 1.79 6.62 11.81 934 1
## sigma 13.34 0.30 12.86 13.81 5213 1par(mfrow = c(1,1))
plot(precis(mTom.2, depth = 2))
# looks a bit worse than the simpler model
# include shelf, using adaptive prior
mTom.3 <- map2stan(
alist(
totleng ~ dnorm(mu, sigma),
mu <- a_species[sp_index] + a_shelf[shelf_index] + b_species[sp_index]*shade,
a_species[sp_index] ~ dnorm(0, 10),
b_species[sp_index] ~ dnorm(0,10),
a_shelf[shelf_index] ~ dnorm(0, sigma_acs),
sigma_acs ~ dnorm(0, 10),
sigma ~ dnorm(0, 10)), data=tom.trim, iter=5000, warmup=2000, chains =2)## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
## Elapsed Time: 7.10618 seconds (Warm-up)
## 5.47194 seconds (Sampling)
## 12.5781 seconds (Total)## The following numerical problems occured the indicated number of times after warmup on chain 1## count
## Exception thrown at line 27: normal_log: Scale parameter is 0, but must be > 0! 4## When a numerical problem occurs, the Hamiltonian proposal gets rejected.## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected## If the number in the 'count' column is small, do not ask about this message on stan-users.##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
## Elapsed Time: 6.25795 seconds (Warm-up)
## 7.16284 seconds (Sampling)
## 13.4208 seconds (Total)## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000207 seconds (Sampling)
## 0.000211 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]plot(mTom.3)## Waiting to draw page 2 of 2
precis(mTom.3, depth=2) # again, n-eff not great## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_species[1] 2.93 4.48 -4.15 10.19 1043 1
## a_species[2] 9.76 4.47 2.74 17.16 1062 1
## a_species[3] 10.36 4.48 3.38 17.58 1053 1
## a_species[4] -5.56 4.53 -12.73 1.69 1124 1
## a_species[5] 14.16 4.49 6.91 21.15 1077 1
## b_species[1] 8.21 4.64 0.94 15.45 1501 1
## b_species[2] -2.67 4.65 -10.12 4.68 1478 1
## b_species[3] 1.25 4.65 -6.36 8.28 1562 1
## b_species[4] 7.06 4.76 -0.50 14.77 1592 1
## b_species[5] 4.75 4.65 -2.89 11.81 1526 1
## a_shelf[1] 50.76 6.20 40.92 60.61 877 1
## a_shelf[2] 57.94 6.18 48.20 67.82 881 1
## a_shelf[3] 53.60 6.19 43.52 63.28 890 1
## a_shelf[4] 32.21 4.43 25.10 39.28 1027 1
## a_shelf[5] 32.95 4.47 25.71 40.05 1045 1
## a_shelf[6] 37.36 4.41 30.31 44.40 1037 1
## sigma_acs 30.15 4.70 22.75 37.39 2170 1
## sigma 14.82 0.33 14.29 15.33 3532 1par(mfrow = c(1,1))
plot(precis(mTom.3, depth = 2))
# looks like shlef has a bigger effect than speciesnow include both acs and shelf in the model
mTom.4 <- map2stan(
alist(
totleng ~ dnorm(mu, sigma),
mu <- a_species[sp_index] + a_shelf[shelf_index] + a_acs[acs_index] + b_species[sp_index]*shade,
a_species[sp_index] ~ dnorm(0, 10),
b_species[sp_index] ~ dnorm(0,10),
a_acs[acs_index] ~ dnorm(0, sigma_acs),
a_shelf[shelf_index] ~ dnorm(0, sigma_acs),
sigma_acs ~ dnorm(0, 10),
sigma_acs ~ dnorm(0, 10),
sigma ~ dnorm(0, 10)), data=tom.trim, iter=5000, warmup=2000, chains =2)## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
## Elapsed Time: 11.0538 seconds (Warm-up)
## 15.1305 seconds (Sampling)
## 26.1844 seconds (Total)## The following numerical problems occured the indicated number of times after warmup on chain 1## count
## Exception thrown at line 25: normal_log: Scale parameter is 0, but must be > 0! 2## When a numerical problem occurs, the Hamiltonian proposal gets rejected.## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected## If the number in the 'count' column is small, do not ask about this message on stan-users.##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
## Elapsed Time: 11.2639 seconds (Warm-up)
## 17.1731 seconds (Sampling)
## 28.4371 seconds (Total)## The following numerical problems occured the indicated number of times after warmup on chain 2## count
## Exception thrown at line 25: normal_log: Scale parameter is 0, but must be > 0! 2## When a numerical problem occurs, the Hamiltonian proposal gets rejected.## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected## If the number in the 'count' column is small, do not ask about this message on stan-users.## Warning in FUN(X[[i]], ...): data with name shelf is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name acs is not numeric and not used## Warning in FUN(X[[i]], ...): data with name species is not numeric and not
## used##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.000225 seconds (Sampling)
## 0.000228 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]plot(mTom.4) ## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
precis(mTom.4, depth=2) # n-eff better now## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_species[1] 15.01 5.92 5.41 24.13 2072 1
## a_species[2] 20.45 6.20 10.69 30.66 2132 1
## a_species[3] 20.98 6.14 11.09 30.52 1935 1
## a_species[4] 8.17 5.83 -1.08 17.53 2032 1
## a_species[5] 23.84 6.26 13.71 33.39 1913 1
## b_species[1] 15.38 4.55 7.76 22.17 2472 1
## b_species[2] 4.68 4.59 -2.52 12.14 2302 1
## b_species[3] 8.31 4.54 0.79 15.22 2404 1
## b_species[4] 15.21 4.73 8.00 23.01 2574 1
## b_species[5] 12.03 4.57 4.92 19.37 2483 1
## a_acs[1] -2.61 5.33 -10.42 6.60 2990 1
## a_acs[2] -7.15 5.05 -15.07 0.95 2617 1
## a_acs[3] 1.29 5.10 -6.54 9.48 2889 1
## a_acs[4] 8.39 5.06 0.13 16.26 2726 1
## a_acs[5] 19.46 5.12 11.08 27.45 2695 1
## a_acs[6] 7.64 5.11 -0.49 15.51 2762 1
## a_acs[7] 15.79 5.07 7.46 23.59 2369 1
## a_acs[8] 11.77 4.84 3.97 19.14 2743 1
## a_acs[9] 8.30 4.82 0.83 15.79 2631 1
## a_acs[10] 6.10 5.19 -2.48 14.06 3349 1
## a_acs[11] 8.13 4.98 0.38 16.31 2378 1
## a_acs[12] -6.64 4.95 -14.66 1.06 2824 1
## a_acs[13] -1.53 4.98 -9.88 6.03 3194 1
## a_acs[14] 3.41 5.01 -4.76 11.06 2738 1
## a_acs[15] -3.95 6.46 -14.26 6.09 3994 1
## a_acs[16] 7.30 4.97 -0.39 15.51 2491 1
## a_acs[17] 9.30 4.79 1.80 16.90 2924 1
## a_acs[18] -9.38 4.88 -16.56 -1.09 2993 1
## a_acs[19] 21.28 5.07 12.88 29.02 2373 1
## a_acs[20] -0.87 4.80 -8.29 6.75 2628 1
## a_acs[21] 6.19 4.83 -1.37 13.88 2671 1
## a_acs[22] 7.83 5.91 -1.84 16.77 3677 1
## a_acs[23] 0.52 5.07 -7.81 7.99 3286 1
## a_acs[24] 0.56 5.12 -7.41 8.75 2885 1
## a_acs[25] 0.05 5.02 -7.60 8.18 2799 1
## a_acs[26] 3.64 5.08 -4.39 11.71 2732 1
## a_acs[27] 6.24 4.80 -1.26 13.99 2868 1
## a_acs[28] -1.18 4.83 -8.68 6.58 2848 1
## a_acs[29] 5.27 4.91 -2.86 12.62 2776 1
## a_acs[30] 11.05 4.84 3.81 19.31 2697 1
## a_acs[31] 15.93 4.97 8.01 23.62 2774 1
## a_acs[32] 4.49 4.97 -3.43 12.44 3175 1
## a_acs[33] 0.17 4.94 -7.89 7.76 3098 1
## a_acs[34] -4.03 4.97 -11.87 3.84 2565 1
## a_acs[35] -1.33 5.13 -8.88 7.27 2284 1
## a_acs[36] -4.60 5.58 -13.60 4.18 3471 1
## a_shelf[1] 26.82 6.35 16.60 36.78 1620 1
## a_shelf[2] 34.70 6.36 25.08 45.28 1613 1
## a_shelf[3] 30.37 6.33 20.27 40.36 1593 1
## a_shelf[4] 15.80 4.50 8.67 22.97 1891 1
## a_shelf[5] 17.34 4.55 9.99 24.54 1913 1
## a_shelf[6] 21.75 4.48 14.45 28.82 1867 1
## sigma_acs 13.06 2.09 9.87 16.38 1918 1
## sigma 13.04 0.29 12.60 13.50 5486 1par(mfrow = c(1,1))
plot(precis(mTom.4, depth = 2))
# compare these 4 models
compare(mTom.1, mTom.2, mTom.3, mTom.4)## WAIC pWAIC dWAIC weight SE dSE
## mTom.4 8083.7 44.3 0.0 1 56.24 NA
## mTom.2 8124.1 39.2 40.4 0 55.24 14.20
## mTom.3 8310.2 14.5 226.5 0 54.42 31.29
## mTom.1 8334.7 10.2 251.0 0 53.84 32.95As I suspected based on the n_eff values, model 4 is much preferred.
Problem 2: Smoking deaths among doctors
In 1961 Doll and Hill sent out a questionnaire to all men on the British Medical Register inquiring about their smoking habits. Almost 70% of such men replied. Death certificates were obtained for medical practitioners and causes of death were assigned on the basis of these certificates. The breslow data set contains the person-years of observations and deaths from coronary artery disease accumulated during the first ten years of the study.
Analyse this data set to determine the posterior probability that smoking increases death by coronary artery disease, that age increases death by coronary artery disease, and that there is an interaction between age and smoking.
You can load the data set and learn about the columns using the commands below
data(“breslow”,package = “boot”) help(“breslow”,package =“boot”) breslow You can think of “person-years” as the number of observations
Note: You almost certainly have the boot package on your computer, but if you do not have the boot package on your computer then you will need to install.packages(“boot”)
Note: do NOT do library(boot). This will make the logit function from boot over-ride the one from rethinking Note: You will probably need to do:
breslow\(n <- as.integer(breslow\)n)
before you can analyze the problem 2 data
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.
Problem 3: Cane Sugar
This data comes from an experiment to measure disease resistance in different varieties of sugar cane.
Is there evidence of differences in disease resistance in the different varieties? Does including an adaptive prior for Block improve the model fit?
You can get the data and learn about it with:
data(“cane”,package=“boot”) help(“cane”,package=“boot”) head(cane) summary(cane)